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CALCULATION OF TEMPERATURE DISTRIBUTIONS IN THIN SHELLS 
OF REVOLUTION BY THE FINITE -ELEMENT METHOD 

By Howard M. Adelman and Donnell S. Catherines 
Langley Research Center 

SUMMARY 

This report describes a procedure for computing the steady-state temperature dis- 
tributions in thin shells of revolution having general meridional curvatures. The pro- 
cedure is based on the finite-element method using a geometrically exact element. 

The procedure is applicable to shells subjected to a very general set of thermal 
loading conditions including applied heat flux at the surfaces, convective heat transfer 
with the surroundings, and equivalent loadings due to specified temperatures on the shell. 
Two general-purpose computer programs based on the procedure have been developed 
and are described. 

The computer programs have been demonstrated for a variety of shell configura- 
tions for which exact solutions were available. These configurations include cylinders 
under a number of different heating conditions, a conical frustum, a truncated hemisphere, 
and an annular plate. Results from the computer programs are compared with the exact 
analytical solutions. 


INTRODUCTION 

Current aerospace applications require the structural analyst to have a knowledge 
of the dynamic response of structural components in a thermal environment. Among the 
most important of these components are thin shells of revolution. As major components 
in aeronautical and aerospace vehicles, shells are subjected to thermal inputs such as 
aerodynamic heating, heating due to contact with an adjacent high -temperature environ- 
ment, and nonuniform heating such as occurs when a small section of the shell is in 
contact with an adjacent heated component. 

A practical method for computing forced response of shells in a thermal environ- 
ment is being developed at Langley Research Center. A generally accepted systematic 
procedure for this problem (outlined below) is being followed: 

(1) Calculate the temperature distribution in the shell due to the thermal inputs. 

(2) Compute the thermal stresses in the shell due to this temperature distribution. 



(3) Compute the natural frequencies and mode shapes of the shell in the presence 
of the thermal stresses. 

(4) Calculate the dynamic response of the shell due to applied mechanical loads by 
modal superposition. 

The first and most basic step in the procedure requires a method for computing 
nonaxisymmetric steady-state temperature distributions in thin shells of revolution 
having general meridional curvature. The resulting temperature distributions would 
then be used as input to other computer programs which would proceed through steps (2) 
to (4). The method used to predict temperature distributions should be fast and accurate, 
and should be applicable to a wide variety of shell geometries and shells having nonhomo- 
geneous properties. 

Two methods have been most successful for solving problems of this type: the 
finite-difference method (ref. 1), and the finite-element method (refs. 2 to 6). The 
finite- element method has the following advantages relative to the finite-difference 
method : 

(1) The finite-difference method operates on the governing differential equation 
whereas the finite-element method operates on an appropriate variational expression 
that is usually somewhat easier to obtain. 

(2) The finite-element method results in a symmetric matrix formulation whereas 
the matrices in a finite-difference formulation are not always symmetric. Symmetric 
matrix formulations are desirable because of the body of mathematical tools available 
for solving the resulting equations. Since computer programs to meet the aforementioned 
requirements did not exist, a new finite-element procedure and computer program were 
developed and are described in this paper. 

A geometrically exact element was employed because of previous success in 
applications to vibration analysis. (See ref. 7.) The computer programs have the capa- 
bility of computing the nonaxisymmetric, steady-state temperature distribution in thin 
shells of revolution subject to a general set of thermal loading conditions including heat 
flux and equivalent loads due to temperatures specified on the shell. They have been 
applied to a variety of shells of revolution. Comparison of results with available exact 
solutions was made and is presented in this paper. 

SYMBOLS 

A matrix which relates temperatures and derivatives of temperature at ends of 

element to coefficients of polynomial representing temperature in element 
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a o>aj,a 2 ,ag coefficients of polynomial representing temperature distribution in 
element 

C column matrix used to specify position of a point where temperature is 

specified and magnitude of specified temperatures (eq. (50)) 

E expression which, when minimized, gives governing equations for steady- 

state temperature distribution in thin shells of revolution 

Fi,F 2 partitions of thermal force vector corresponding to unknown and known 
temperatures, respectively 

f thermal force vector 

G expression which, when minimized, gives governing equations for tempera- 

ture distribution for shell in which temperatures are specified at discrete 
points 

H convective heat-transfer coefficient 

h thickness of shell 

K number of finite elements used to represent shell 

k(s,z) thermal conductivity 

r h/2 

K(s) = \ k(s,z)dz 

J -h/2 

k element number 

L length of shell 

M number of circumferences on which temperatures are specified 

N largest Fourier index in a set of calculations 

n Fourier index 

P number of discrete points at which temperature is specified 
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Q column matrix whose elements are coefficients of linear terms in E 

q heat flux 

R matrix whose elements are coefficients of a quadratic form in E (eq. (5)) 

r perpendicular distance from shell axis to a point on shell middle surface 

S conductivity matrix 

S 11> S 12> S 21> S 22 partitions of S (eq. (40)) 
s meridional coordinate 

T specified temperature 

U vector containing temperatures at element junctures 

u temperature 

V column matrix used to specify position of a point where temperature is 

prescribed and also the magnitude of the temperature (eq. (51)) 

X matrix which describes polynomial approximation of temperature dis- 

tribution in an element (eq. (14)) 

x meridional coordinate in an element 

y vector defined in equation (12) 

z coordinate normal to shell middle surface measured from middle surface 

column matrix whose elements are the coefficients of the polynomial approxi- 
mating the temperature distribution in an element 

(k # K) 

(k = K) 

(k* 1) 

(k = 1) 
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e 


meridional length of an element 


6 circumferential coordinate 

X Lagrange multiplier 

Subscripts: 

F first edge of shell 

I inside surface of shell 

k element number or juncture number 

L last edge of shell 

m index referring to circumference at which temperature is specified 

O outside surface of shell 

p index which refers to point at which temperature is specified 

cn Fourier cosine component 

sn Fourier sine component 

Superscripts: 

cn Fourier cosine component 

sn Fourier sine component 

T transformed matrix 

Primes indicate differentiation with respect to x or s. 

ANALYSIS METHOD 

In this section of the report, the governing equations are derived. The following 
procedure is followed: A functional is obtained which, when minimized, yields the 



boundary-value problem governing steady-state conduction of heat in a thin shell of 
revolution. This functional is obtained by specializing a more general expression 
(refs. 2 to 6) by making the following assumptions: 

(1) There is no time variation of temperature. 

(2) Temperature gradients through the thickness of the shell may be ignored. 

(3) The geometry and properties of the shell are independent of the angular 
coordinate 6. 

The functional constitutes the basis for the entire approach. Once it is derived, the 
analysis follows a well-known procedure described in references 2 to 6. 

Each variable in the functional is expanded in a Fourier series around the shell 
circumference so that the calculation of temperature distribution reduces to computing 
only the meridional variation of temperature. In accordance with the finite-element 
method, the shell is represented by a number of elements in which the temperature is 
approximated by polynomials. In the present procedure the elements are geometrically 
exact in that no approximation of the shape of the shell is necessary. Conductivity 
matrices and thermal force vectors are derived for each element and these element 
matrices are superimposed to form the corresponding matrices for the complete shell. 

It is assumed in this analysis that both the temperature and its first derivative are con- 
tinuous at junctures between adjacent elements. According to reference 6, as a mini- 
mum, only temperature need be continuous in order to guarantee convergence of the 
method. Nonetheless, specifying continuity of the first derivative is a reasonable condi- 
tion to impose since in most of the anticipated applications of the present method, the 
temperature distributions will be smooth functions. In the cases where the actual tem- 
perature distributions have slope discontinuities, the method based on continuous slopes 
will still converge since the conditions expressed in reference 6 are satisfied. The 
resulting solution will cause the "corner” at the discontinuity to be rounded. 

The matrix equations for determining the temperature distribution are derived by 
minimizing the approximated functional with respect to the generalized coordinates which 
in this analysis are the values of temperature and first derivative of temperature at the 
ends of each element. 

For the purposes of the following development, reference is made to figure 1. In 
this figure u(s,0) represents the temperature of a point on the shell defined by the 
meridional coordinate s and the circumferential coordinate 6. The ambient tempera- 
tures of the media outside of the first and last edges of the shell are represented by 
Up(0) and u l (6»), respectively. The ambient temperature of the medium at the outer 
surface of the shell is Uq(s,0) and the ambient temperature of the medium at the inner 
surface of the shell is Uj(s,0). Heat flux is applied to the outer and inner surfaces of 
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q F 0>) 

Figure 1.- Nomenclature associated vith transfer of heat 
in a thin shell of revolution. 

the shell by sources qQ(s,0) and qj(s,0), respectively. The applied heat flux at the 
first and last edges are given by q-p(0) and q^#), respectively. The perpendicular 
distance from the axis of revolution to a point on the shell middle surface is given by 
r(s). The distance from the origin of the meridional coordinate s to the first edge is 
s Q and the distance to the last edge is L. The thickness of the shell is given by h(s). 

The functional for the temperature distribution for bodies of general shape is well 
known. (See, for example, refs. 2 to 6.) Specialization of the functional for steady-state 
heat conduction in isotropic thin shells of revolution gives the following functional: 


E “ 1 1ST k(s ’ z) |(S) + r(s) ds de dz + ii* H o< s >(V - u o u ) r(s) ds de + Jj* h i< s >(y ' u i u ) r(s) ds de 


q o( s >^ 


ur(s) ds de - j j qj(s,0) ur(s) ds de 


q F (e) ur(s 0 ) dz de 


ur(L) dz de 
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where 


(1) k(s,z) is the thermal conductivity of the shell material 


(2) Hq(s) and Hj(s) are the convective heat -transfer coefficients at the outer 

and inner surfaces of the shell 

(3) Hp* and are the convective heat -transfer coefficients at the first and 

last edges of the shell. The integrations indicated by \ \ and \ \ 

JjQ-p JJ( 

are taken over the first and last edges, respectively, of the shell. 

Radiation effects are entirely neglected in this analysis. 

For a shell of revolution, the variables may be represented in terms of their 
Fourier components as follows: 

■'n 


CL 


n=0 
N 


N N 

u(s,0) = ^ u cn (s) cos ne + ^ u sn (s) sin ne 

n=l 


n=0 


N 

u O (s ’ 6 > = Z U 0,cn^) cos ne + Z U 0,sn^ s ) sin n6 

n=l 


N N 

U I< S , 0 )= £ u I>cn (s) cos n0 + £ u I>gn (s) sin n 6 
n=0 ’ n=l 


N N 

q o ( s ,0) = £ \),cn (s ) cos n6 + I <lo, S n (s ) sin n6 

n=0 n=l 

N N 

q i(s,e) = ^ qj cn (s) cos ne + ^ q Isn (s) sin *6 
ri=0 ’ ri=l 


N N 

q F (e) = £ q FiCn cos + £ q Fsn sin nf> 
n=0 n=l 

\ 

(Equations continued on next page) 
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cos nfi + 


( 2 ) 


N 

q L (e> - £ q L . 


_ , cn Li L.sn 

n=0 n=l 


N 

I q L,£ 


sin n0 


N IN 

n F (e) = ^ u F>cn cos ne + £ u F>gn sin n0 
n=0 


N 

l 


N IN 

U L^> = £ U L,cn cos ne + £ u L,sn sin n0 
n=0 


N 

i 


Substitution of equations (2) into equation (1) and integration with respect to 6 
from 0 to 2tt yields the following form for the functional E: 

N N 



n=0 n=l 


The forms of E cn and E sn are identical. The formulation of the finite-element 
representation of E sn is the same as the formulation for E cn . Therefore, in order to 
avoid needless repetition, the derivation will be carried out for E cn with the under- 
standing that the treatment of E sn is the same. 

The form of E cn is 


r* ph(s)/2 

Ec " = fW- h (s)/2 k(S ’ Z) 


( U cn) 2 + Jl% rdz ds + 7r J - u 0;Cn u cn j r ds + * J - u IjCn u cn ]r ds 


r* r » r* h F/^ l / 2 

- * J q O,cn u cn r ds * * J 1l,cn u cn r ds ' n J. h ^ 2 lF,cn u cn r ( s o) dz - n J ^.cn^n^L) dz 

r h F /2 / u 2 \ ^hr/2 /u 2 \ 

+ * J-h F /2 Hf V“T ' u cn u F,cnJ r ( s o> d * + * /2 H L^ ' W L ,cn) r < L > 


dz 


(4) 


For n = 0 the value of E cn is twice the value indicated by equation (4). The integra- 
tion denoted by ^ is taken over the length of the shell meridian and h F and hL are 
the shell thicknesses at the first and last edge of the shell, respectively. 
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Representation of Shell by Finite Elements 

The shell is represented by geometrically exact elements by which it is meant that 
the element shape coincides with the part of the shell being represented by the element. 
A typical idealization is shown in figure 2. By counting elements from the first edge of 
the shell, the following definitions are made: 

K total number of elements 

e k length of kth element measured along meridian of shell 

x meridional coordinate inside kth element, measured along meridian from 

- 6if ^ k 

center of kth element so that ^ x % — 

2 2 

s k distance along meridian from first edge of shell to center of kth element 

s 0 distance from origin of s to first edge of shell 

From the foregoing definitions, it follows that 
s = s 0 + s k + x 

Quantities evaluated at the first edge of the kth element (that is, x = -e k /2^ are 
denoted by the subscript k. Quantities evaluated at the second edge of the kth element 
^that is, x = e^/2^ are denoted by the subscript k + 1. 



Figure 2.- Typical idealization of shell of revolution 
by geometrically exact finite elements. 
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When discussing the complete shell as idealized by finite elements, it is necessary 
to refer to the term "element junctures." The kth element juncture is defined to be the 
circumference which is the first edge of the kth element. The one exception to this 
general definition is that the last edge of the shell is designated as the K + 1 element 
juncture. (See fig. 2.) 

By following the conventional steps of the finite-element procedure, the integrations 
are taken to extend over a single element. 

After integration with respect to z, equation (4) may be written in the following 

form: 



( 7rh F r F H F u F>cn ^u l cn - (^h L q L}Cn r L ^u K+1>cn 


" |( irh L H L r L) u K+l,cn + ( lh L H L r L u L,cn) u K+l,cnJ 6 K 

where ^R^ cn |j is a two-by-two symmetric matrix where elements are functions of x 
as given below: 


(5) 


["k” n) ] - 


2 

+ h g (x) r(x) + Hj(x) r(x) 0 


r 2 (x) 


K(x) r(x) 


( 6 ) 


r h/2 

K(x) = \ k(x,z) dz 
J -h/2 


(7) 


«cn <x > = , So,cn< x )--%,cnW 


( 8 ) 
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r-p = r(s D ) sr at first edge of shell 
r L = r(L) = r at last edge of shell 
U 1 cn tem P erature component at first edge of shell 
U K+1 cn tem P eratur e component at last edge of shell 



(k* 1) 
(k = l) 



(k * K) 
(k = K) 


As an approximation, the temperature is assumed to have the following polynomial 
form within the kth element: 


u cn M = af? ♦ a< c ">x + a< c ”>x 2 + a< c f)x 3 
tn U > K l,k 2,k 3,k 

Equation (9) may be written as 

, T 


where 


u cn (x) = {,('»} (y} 


( y (cn)} = ; 


r a ( Cn )^ 

a 0,k 

a ( cn ) 

1, k 

a M 

2, k 

a ( cn ) 

v. 3 > k .> 






X 

c 

X* 


( 9 ) 


( 10 ) 


(ID 


( 12 ) 
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The variables sS^ n \ a[ c , n \ a^°P\ and a^ c ^ are the undetermined coefficients asso- 
0,k ’ l,k 2,k 3,k 

ciated with the nth Fourier cosine component in the kth element. 

The relationship between the temperature and its first derivative at the ends of an 
element and the coefficients of the assumed polynomial is derived as follows. 

Using equation (9) and a differentiation of both sides of the equation yields 



(13) 


(14) 


Inserting x = -e^/2 and x = e k/^ in* 0 appropriate locations in equation (13) and 
inverting the resulting matrix equation gives the following relation: 


where 


(y (cn) ) = [ A k]{u 




6 k/ 8 

1/2 


-1/4 

3/2e k 

-1/4 

V 2e k 

0 

V 2£ k 

i/2 

-V4 



(15) 


(16) 



(17) 
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in terms of the 


Substituting equations (10) and (13) into equation (5) gives E cn k 
undetermined coefficients of the polynomials as 

E cn,k * | M ^ (Cn> > * 


Define 



bp - irhprpHp 

B L n) ’ ’(Vl-Kcn + h L“L,c„)] 
b L = 7rh L r L H F 


(18) 

(19) 

( 20 ) 

(21) 

( 22 ) 

(23) 

(24) 
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By utilizing equations (19) to (24) , equation (18) may be written as 


E c „,k - Kr(“)> T [z‘“>]{r(-)> - {r< cn >} T (<£">} 

- [bf”' u l,cn ■ | b F"l,cn] 5 l ■ (®L C ”' u K+l,cn " | b L u K+l,cn] 5 K * 25) 
Substituting equation (15) into equation (25) yields 

E c„, k - i(4 cn j} T [ A k] T [4 cn) ] m H cn> } - H cn D T M T Hi 

- Uj cn " | b F“l, C n] 6 l ' u K+l,cn " | b L u K+l ,cn] 6 K < 26 > 


The conductivity matrix and thermal force vector may be identified for each ele- 


ment as 


4 cn >] and {4 C ”>}, 


respectively, where 


4 cn) ] - |>] T [4 cn> ] M ♦ [ b 0^ ♦ [W]4 


(27) 




5 k 


(28) 


where 



0 0 0 

0 0 0 
0 0 0 


0 0 0 0 


(29) 
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0 

0 

0 

0 


0 

0 

5 l 

0 


0 

0 

0 

0 



( 30 ) 


(31) 


(32) 


In view of the numbering convention adopted for the elements, the second edge of 
the kth element coincides with the first edge of the k + 1 element. In this analysis it 
is assumed that the following conditions of compatibility hold at each such juncture: 


f. 


(cn) 

k+1 


U M’ 

l k+1 ^ 




''(cup 

k+1 


u 


k 

element 


U M' 

v > +1 J 


k+1 

element 


(1 £ k ^ K) (33) 


These compatibility conditions require that the temperature be continuous with a continu- 
ous first derivative at all points along the shell. The equality between the derivatives is 
strictly valid only for cases where the heat flux and thermal conductivity are continuous. 

The functional E cn associated with the nth harmonic for the entire shell is given 
by the following summation: 
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( 34 ) 


K 

E cn = Y E cn,k 
k=l 


where E cn k is given by equation (26). 

Use of equation (33) in equation (34) gives 


E cn = - {u (cn) } T [s (cn )] {u< cn >} - {u (cn )} T {f ( cn)} (35) 

2 


where 

s( cn ) thermal conductivity matrix which is a symmetric positive definite 

matrix of order 2(K + 1) x 2(K + 1) (see ref. 2) 

f( cn ) thermal force vector which is of order 2(K + 1) x 1 

u( cn ) a vector containing all temperature components and derivatives 

Thus, 


(cn) 




{u< cn >} = I 


U 


Tj( Cn ) 

u 2 


u( cn ) 

v. U K J 


2(K+l)xl 


(36) 


In the present procedure, the matrices S^ cn ^ and f( cn ) are constructed by the 
well-known procedure of superposition of element matrices as illustrated in figure 3. 
The superposition to form S^ cn ^ consists of placing the element matrices together 
so that the matrix elements in the lower 2x2 block of the kth matrix add to the corre- 
sponding matrix elements in the upper 2x2 block of the k + 1 matrix. In the case of 
f( cn ), the superposition consists of placing the element vectors together so that the last 
two numbers in the kth vector add to the first two numbers in the k + 1 vector. 
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Figure 3.- Illustration of superposition of element matrices. 


By recalling that the treatment of E sn is identical to the treatment of Ecn> 
equation (35) is substituted into equation (3) and the following expression is obtained: 

E = < Cu^ cn )} >T [s( n )J - (u( cn )} T {f( cn )}) 

n=0 

+ £ (|frW> T [s( n )]{u( sn )) - (uM> T { f( sn )}j 

n=l x 

Since ^ s (cn)J 

is identical to jj3^ sn )J , they are both denoted by [s^ n ^j . 


Formulation of Governing Matrix Equations 

The equations for determining the vectors U^ cn ^ and tK sn ) are derived by 
minimizing E with respect to the variables in these vectors. It is seen from equa- 
tion (37) that there are no cross products between sine and cosine vectors. Further, 
there are no cross products between different sine components or between different 
cosine components. Therefore, the minimization of E with respect to a given sine or 
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cosine component gives an equation which involves only that component. Thus, each 
cosine component and each sine component may be calculated independently. This prop- 
erty is of practical significance and its exploitation will be illustrated in the discussion 
of the computational method. In the following derivations, the operations for the typical 
cosine component are performed.with the understanding that all the cosine components as 
well as all the sine components are treated in an identical fashion. 


Formulation for case where temperature over entire shell is unknown .- If the tem- 
perature at every point on the shell is unknown, the equations for determining a given 
Fourier component of temperature are found by minimizing E with respect to each of 
the elements in the vector U^ cn \ The minimization is equivalent to the following set of 
equations, 2(K + 1) in number: 


3E 


9u 


(cn) 


= 0 


8E 


3u 


(cn)’ 


= 0 


(k = 1, 2, . . ., K + 1) 

> 

(k = 1, 2, . . ., K + 1) 


(38) 


Equations (38) can be expressed in the form: 

[sW] {u( cn) } = {f^} (39) 

Equation (39) may be solved for {u( cn )} by inverting the matrix since it is 

positive definite and therefore nonsingular (ref. 2). 

Recall that the elements in the vector {u^ cn ^ are values of the Fourier compo- 
nents of temperature and meridional derivatives thereof at each element juncture. By 
using equation (15), the coefficients of the polynomial within each element may be com- 
puted from the elements in {\j( cn )}>. Then the temperature component within each ele- 
ment may be computed at as many points as desired by use of equation (9). 


Formulation for case where temperature is specified on circumferences .- There 
is no loss in generality in assuming that the circumferences along which temperature is 
specified are element junctures. When this assumption is made, the vector U^ cn ^ may 
be partitioned into an unknown part v^ cn ^ and known part w^ cn ' . Then 


E 


cn 



? (cn) 

’ll 


21 


o(cn) 

b 12 


s< c . n > I s< cn) 


22 




(40) 
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where 


S M a square matrix whose rows and columns correspond to the unknown elements 

of U^ Cn ^. It is of order 2(K + 1) - M 

S M a rectangular matrix having 2(K + 1) - M rows and M columns. The 

rows correspond to unknown elements of U^ cn ^ and the columns corre- 
spond to known elements of l/ cn \ 

transpose of 

sg» a square matrix of order M whose rows and columns correspond to known 

elements of u^ cn ^ 

f( C 11 ^ a partition of f( cn ) whose rows correspond to the unknown elements 

1 ot U (cn > 

F^ cn ^ a partition of f^ cn ) whose rows and columns correspond to the known 

elements of U^ cn ^ 


If the temperature of the shell is specified on the edges of the shell, then the variables 
U 1 cn anc * u k+l cn are deluded i n the vector w^ cn \ The treatment of edge conditions 
is discussed in more detail in a later section. Minimizing E with respect to the ele- 


ments in v 


(cn) 


yields the equations for determining the unknown temperatures. The 


minimization is the same as that of equation (38) except the differentiations are taken 
only with respect to the elements of v^ cn \ The resulting matrix equation may then be 
written as 


[sjj n, ]{!v< cn >} = - [s<->] {w< cn >} 

This equation is solved by inverting the matrix j"s^ n J . This matrix is nonsingular 
because it is a submatrix of the positive definite matrix 


Formulation of equations for case where temperatures are prescribed at discrete 
points .- For the purpose of the following development, the reader is referred to figure 4. 
The pth point at which temperature is specified is located by the coordinates Sp,0p. The 
temperature at that point is denoted Tp. The number of the element in which the pth 
point is located is denoted by k p and the corresponding value of given by sjqp. 

The meridional distance from the center of the kpth element to the pth point is denoted 
by x p . 
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Figure 4.- Nomenclature associated with specification of temperature 
at a discrete point on a shell. 

From these definitions, the following relations are obtained: 


u ( s p>^p) ~ T p 

(P = 1, 2, . 

• P) 

(42) 

X P = s p ■ s kp " s ° 

(P= 1, 2, . 

• P) 

(43) 


where P is the number of discrete points at which temperature is specified. 

The derivation of the governing equations for this class of problems is facilitated 
by the use of Lagrange multipliers. It is required to find stationary values of E (eq. (3)) 
subject to the constants of equation (42). 

It is therefore necessary to find stationary values of the functional G where 

G = E + - T X J + . . . + Ap|u(sp,e p ) - Tpj (44) 

and E is given by equation (37) . The quantities Aj, . . A p are Lagrange multi- 
pliers. Using the first of equations (2) with equations (10) and (42) gives the following 
relation: 

N T N T 

U ( S P’ e p) = ^{r* Cn ?} {y p } cos n0 p + ^ {y( sn )} {y^ Sin n0p (45) 

n=0 n=l 
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where 


r l^ 


C y p} - { 


4 

x p 

V J 


Substituting equation (15) into equation (45) gives 

-(Vp) - 1 ( u fe n) ) T {v + 1 (<"5 Ot p > 

n=0 n=l v 


where 

(V> = [ A kp] T (yp} cos ne p 

T 

{*kp} =[ A kp] {y p } sm nS p 


Two vectors 



are defined as follows: 


fol 


T 

2(kp - 1) rows 




L°J 


fol 


2(kp - 1) rows 



= a kp )_i_rows 


L°J 


(46) 


(47) 


(48) 

(49) 


(50) 


(51) 
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I 


Equation (47) may be written with the use of equations (48) to (51) as 

N 

u(s p ,0p) = £ (u< Cn )} T {cg) + ^uM>(vg} (p = 1, . . P) (52) 

n=0 ^ J ^ J 

Use of equation (52) in equation (44) gives required form of G: 

G = E + X 1 Y {u (cn) } T f cjjjl + \'} j t u(Sn) ^ T (vQ| - X 1 T 1 
n=0 ^ ^ n=l ^ .< 


+ Xt 


, | o (u (cn) > T {cg} 


+ X-, 


2(^")} T {vW} - 


Xp T p 


(53) 


The governing matrix equations for this problem are derived by minimizing G 
with respect to the elements in U^ cn ^ and U^ sn ^ and with respect to the Lagrange 
multipliers to Xp. In contrast to the problem considered previously, all the har- 
monics must be computed simultaneously. The minimization of G implies the following 
conditions: 


8G =(T 


da 


(cn) 

k 

3G 


da 


(cn) 

k 

3G 


- = 0 


da 


(sn) 

k 

3G 


= 0 


^k 

3G 

3X„ 


(sn)’ 


= 0 


= 0 


(n = 0, 1, . . N) 
(k = 1, 2, . . K + 1) 


(n = 1, 2, . . . , N) 
(k= 1, 2, . . K + 1) 


(p = 1, 2, . . P) 


(54) 


The number of equations represented by equations (54) is 4(K + 1)(N + 1/2) + P. 
These equations may be written as follows: 
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[s^{u(e n )> - (f(cn)} +X 1 {cgj + . . . +X p |^C^J =0 

(n = 0, 1, . . N) (55a) 

[s (n |] (u( sn )} - {f( sn )} + Xi(vg) + . . . + {^ v kp| =0 

(n= 1, 2, . . N) (55b) 



(p = 1, 2, . . P) (55c) 

The actual form of this set of equations written as a single matrix equation is 
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The matrix in equation (56) is nonsingular (ref. 2) so that the solution of the equa- 
tion may be effected by inverting the matrix. The resulting vector can be broken down 
into the subvectors U^ cn ^ and U^ sn ^ for each value of n. Coefficients of the poly- 
nomials may be determined by use of equation (13) and the polynomial representations 
of U cn (x) and U gn (x) within each element may be computed. 

If the temperatures, in addition to being specified at discrete points, are also speci- 
fied along circumferences, the only changes in equation (55) are as follows: 

(s^ n |] is replaced by j^S^j 

(f( cn )} is replaced by - [sjg]' {w( cn )} 

(f( sn )} is replaced by {w^ sn )} 

<(jj( cn f)> i s replaced by (jr( cn fy 
{u (sn) } is replaced by {v^ sn ^} 

OUTLINE OF COMPUTER PROGRAMS 

In this section of the paper, the basic computing procedure is illustrated by means 
of flow charts. In addition, the input to the computer programs is described along with 
the form of the output. The section is divided into two parts. The first part deals with 
the computer program which treats the case where temperatures, when specified on the 
shell, are specified only on circumferences and the second deals with the program for 
the case where temperatures are specified at discrete points as well as on circumferences. 

Computing Method for Calculations Wherein Known Temperatures 
are Specified on Complete Circumferences Only 

The organization of the computing steps for this type of calculation is illustrated 
in block diagram 1. 
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The inputs are of three types: input independent of the circumferential harmonic, 
cosine -type input, and "sine type" input. The input which is independent of the circum- 
ferential harmonic is as follows: 

K number of elements used to represent shell 


e k 


So 


r(s) 

h(s) 


K(s) 


Hrfs) 

Hq(s) 

NBEG 


length of kth element where k - 1, 2, . . ., K 
meridional distance from origin of s to first edge of shell 
shell radius 
shell thickness 

integral of thermal conductivity through shell thickness (see eq. (7)) 
convective heat -transfer coefficient at inner surface of shell 
convective heat -transfer coefficient at outer surface of shell 
initial value of n 


NIjAST final value of n 

M number of element junctures at which temperature is specified 

m(i) an array in which the ith entry is the number of the ith element juncture at 

which temperature is specified, i = 1, 2, . . M 

Hp convective heat -transfer coefficient at first edge of shell 

h l convective heat-transfer coefficient at last edge of shell 

ININ number of integration intervals used to carry out the evaluation of matrices 

Z and Q by the trapezoidal rule 

IPL$T number of locations within each element that temperature distribution 
is printed 
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The integrations required for computing matrices Z and Q are performed 
numerically by using the trapezoidal rule with 100 integration intervals. The tempera- 
ture distributions are computed and printed at 10 locations in each element. Both of 
these numbers may be changed by varying the input quantities ININ and IPL0T. 

The cosine-type input is as follows: 


u O,cn (s) 

W s > 

q F,cn 

u F,cn 


nth cosine harmonic of temperature of medium in contact with inside of shell 

nth cosine harmonic of temperature of medium in contact with outside of shell 

nth cosine harmonic of net heat flux applied to surface of shell 

nth cosine harmonic of heat flux applied at first edge of shell 

nth cosine harmonic of temperature of medium in contact with first 
edge of shell 


cn nth cosine harmonic of heat flux applied at second edge of shell 

u T nth cosine harmonic of temperature of medium in contact with second 

L,cn 

edge of shell 

T m cn nth cosine harmonic of specified temperature at element juncture m(i) 
(where i = 1, 2 , . . ., M) 


The sine -type input is as follows: 


u I,sn< s > 

u O,sn (s > 


<W S > 


q F,sn 

F,sn 

q L,sn 


u 


nth sine harmonic of temperature of medium in contact with inside of shell 
nth sine harmonic of temperature of medium in contact with outside of shell 
nth sine harmonic of net heat flux applied to surface of shell 
nth sine harmonic of heat flux applied at first edge of shell 

nth sine harmonic of temperature of medium in contact with first edge of shell 
nth sine harmonic of heat flux applied at last edge of shell 
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nth sine harmonic of temperature of medium in contact with last edge of shell 


u» 

L,sn 

T m sn nth sine harmonic of specified temperature at element juncture m(i) 
(where i = 1, 2, , . M) 


Edge conditions must be specified at each edge of the shell. These edge conditions 
can be any of the following types: 

(1) specified temperature at edge 

(2) applied heat flux at edge 

(3) convective edge condition 


(4) insulated edge condition 

The appropriate edge condition at an edge is specified by use of the input quantities 


H 


F’ 


H 


L> ^FjCn’ ^F,sn’ ^L,cn’ ^L,sn’ u F,cn’ u F,sn’ u L,cn ; 


and u 


L,sn' 


The 


method of specifying edge conditions is summarized in table I. 


TABLE I.- METHOD OF SPECIFYING EDGE CONDITIONS 


Edge condition 

Method of specifying edge condition 

Specified temperature at first 

tf F = 0 

edge 

q F,cn = q F,sn = u F,cn = u F,sn = 0 
T l,cn and T l,sn specified 

Specified temperature at last 

o 

ii 

W 

edge 

q L,cn = q L,sn “ u L,cn = u L,sn “ 0 
T K+l,cn and T K+l,sn specified 

Applied heat flux ai first edge 

q F,cn and q F,sn 
H F = u F,cn = u F,sn = 0 

Applied heat flux at last edge 

q L,cn and q L,sn s P ecifit ' d 
H L = u L,cn = u L,sn “ 0 

Convective edge condition at 

q F,cn = q F,sn = 0 

first edge 

Hp specified 

u F,cn and ‘ J F,sn s P ecified 

Convective edge condition at 

q L,cn " q L,sn = 3 

last edge 

h L. u L,cn : u L,sn s P eci£ied 

First edge insulated 

q F,cn = q F,sn = 0 

u F,cn “ u F,sn 
H f = 0 



Last edge insulated 

q L,cn ” q L,sn = ® 

u L,cn ’ u L,sn " 0 
H l - 0 

Applied heat flux combined with 

q F,cn’ q F,sn specllied 

convection at first edge 

u F,cn; u F,sn s P eciJiud 
II F specified 

Applied heat flux combined with 

q L,cn; q L,sn s P ecUied 

convection at last edge 

u L,c n ; u L.sn specilied 
H l specified 
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The output of this computer program consists of tabulations of the functions U cn (s) 
and U sn (s) for all values of n in the calculations. 

Computing Method for Calculations Wherein Known Temperatures 
are Specified at Discrete Points as Well as on 
Complete Circumferences 

The organization of the computing steps for this type of calculation is illustrated in 
block diagram 2. 

The input to this program includes all the input of the previous program plus the 
following additional input: 

P number of discrete points at which temperature is specified 

Sp meridional coordinate of pth point where temperature is specified 

(p = l,2, . . . , P) 

0p circumferential coordinate of pth point where temperature is specified 

(p = 1 , 2, . , . , P) 

Tp temperature specified at pth point (p = 1, 2 , . . ., P) 

kp element number in which pth point is located (p = 1, 2, . . ., P) 

The output of this program consists of tabulations of the functions Ucn( s ) an d 
U sn ( s ) f° r the range of n considered. 
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APPLICATION OF METHOD 


i 


Description of Shells Analyzed 

In order to demonstrate the versatility of the present method and the accuracy of 
the results obtained, a number of sample calculations were made. Seven shell configura- 
tions were analyzed. For six of these, exact solutions were used for comparison with 
results of the present method. These solutions were either available from previous 
investigators or are derived in the appendix. The seventh configuration was of interest 
but no exact solution was available. The shell configurations were as follows: 

(1) An annular disk insulated at the outer edge, held at fixed temperature at inner 
edge v/ith convection from the lateral faces. Results are compared with the exact solu- 
tion from reference 8. 

(2) A cylindrical shell v/ith both edges at fixed temperature with convection to a 
medium at a uniform ambient temperature. Results are compared with an exact solu- 
tion derived in the appendix. 

(3) A conical frustum subjected a nonaxisymmetric heat flux over the outer surface. 
Edges are held at fixed temperatures. Results are compared with an exact solution 
derived in the appendix. 

(4) A truncated hemisphere subjected to a uniform heat flux at its outer surface. 
Edges are kept at fixed temperatures. Results are compared with an exact solution 
derived in the appendix. 

(5) A cylindrical shell with the circumference at midspan of the cylinder held at a 
fixed temperature and edges held at a different fixed temperature. Results are compared 
with the exact solution in reference 9. 

(6) A cylindrical shell with a generator held at a fixed temperature. Edges of the 
shell are insulated. Results are compared with the exact solution in reference 9. 

(7) A cylindrical shell with a single point held at fixed temperature. Edges are 
insulated. 

Dimensions of the shells and numerical values of the parameters describing the 
thermal environment are assigned convenient values to simplify the calculations. Param- 
eters describing the shells analyzed are listed in table II. The corresponding input to the 
computer programs are listed in table III. As indicated in table III, configuration 5 was 
analyzed by using two different finite-element representations, one which had 10 equally 
spaced elfements and another which had 20 elements with more elements near the center 
of the cylinder than near the ends. 
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TABLE H.- SAMPLE CALCULATIONS 


Parameter 

Annular 

plate 

Cylinder: 

heated 

edges 

Conical 

frustum 

Truncated 

hemisphere 

Cylinder: 

heated 

midlength 

Cylinder: 

heated 

generator 

Cylinder: 

heated 

point 


n = 0 

n = 0 

tn 

VII 

e 

VII 

o 

n = 0 

n = 0 

0 5 n 5 4 

VII 

c 

VII 

o 

h(s) 

0.001 

0.001 

0.001 

0.001 

0.001 

0.001 

0.001 

K(s) 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

r(s) 

s 

1.0 

)j3/2s 

cos(s) 

1.0 

1.0 

1.0 

r'(s) 

1.0 

0.0 


-sin(s) 

0.0 

0.0 

0.0 

H 0 (s) 

0.5 

1.0 

0.0 

0.0 

1.0 

1.0 

1.0 

H x (s) 

0.5 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

u Q (a,e) 

0.0 

1.0 

0.0 

0.0 

1.0 

0.0 

0.0 

Uj(s,0) 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

q 0 ( s > e > 

0.0 

0.0 

cos n0 

1.0 

0.0 

0.0 

0.0 

qjts,©) 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Condition at first edge 

u = 1 

u = 0 

u = 1 

u = 1 

u = 1 

u’ = 0 

u' = 0 

Condition at last edge 

u’ = 0 

u = 0 

u = 6 

u = 2 

u = 1 

u’ = 0 

u’ = 0 

Additional condition 

None 

None 

None 

None 

u(0.5) = 0 

u = 1 at 0 = 0 

u(0.5, 0) = 1.0 

Meridional length, L 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 


TABLE III.- INPUT TO PROGRAM FOR SAMPLE CALCULATIONS 


Parameter 

Annular 

plate 

Cylinder: 

heated 

edges 

Conical 

frustum 

Truncated 

hemisphere 

Cylinder: 
heated mkllength, 
10 elements 

Cylinder: 
heated midlength, 
20 elements 

Cylinder: 

heated 

generator 

Cylinder: 

heated 

point 

K 

5 

5 

5 

5 

10 

20 

5 

10 

e k 

a 0.1 

a 0.2 

a 9/5 

a C.2 

a 0.1 

(b) 

a 0.2 

a l/6 

s 0 

0.5 

0 

1.0 

0.0 

0.0 

0.0 

0.0 

0.0 

r(s) 

S 

1.0 

(t/V 2)s 

cos(s) 

1.0 

1.0 

1.0 

1.0 

h(s) 

0.001 

0.001 

0.001 

0.001 

0.001 

0.001 

0.001 

0.001 

K(s) 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

H^s) 

0.5 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

H q (s) 

0.5 

1.0 

0.0 

0.0 

1.0 

1.0 

1.0 

1.0 

NBEG 

0 

0 

0 

0 

0 

0 

0 

0 

NLAST 

0 

0 

0 

0 

0 

0 

4 

4 

M 

1 

2 

2 

2 

2 

2 

0 

0 

m(i) 

1 

1, 6 

1, 6 

1, 6 

1, 6, 11 

1, 11, 21 

— 

— 

H F 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

h l 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

u I,cn (s) 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

u O,cn< s > 

0.0 

1.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

q cn (s) 

0.0 

0.0 

1.0 

1.0 

0.0 

0.0 

0.0 

0.0 

q F,cn 

0.0 

0.0 

0.0 

0.0 

0.0 

! 0.0 

' 0.0 

0.0 

u F,cn 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

q L,cn 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

u L,cn 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

T m,cn 

1.0 

0.0, 0.0 

1.0, 6.0 

1.0, 2.0 

1.0, 10.0, 1.0 

1.0, 10.0, 1.0 

0.0 

0.0 

u I,sn< B > 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

u O,sn (s) 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

q sn< s > 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

q F,sn 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

u F,sn 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

q L,sn 

0.0 

0.0 

0.0 

i o.o 

0.0 

0.0 

0.0 

0.0 

u L,sn 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

T m,sn 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

p 







1 1 

1 
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(c) 

A 

II 














<d) 

0.5 

s p 














a 0.0 

0.0 








T n 







a 1.0 

1.0 

1 P 

t 







a All equal. 

b €j to €g = 0.08, e 6 to e 15 = 0.0?., e 16 to € 20 = 0.08. 

c kj = 1 k 2 = 1 k 3 = 2 k 4 = 2 k 5 = 3 kg = 3 k 7 - 4 kg » 4 kg * 5 k 10 =» 5 k n =. 5. 

d s 1 = 0 & 2 = 0.1 s 3 = 0.2 s 4 = 0.3 s 5 = 0.4 Sg = 0.5 s ? = 0.6 Sg = 0.7 s fl = 0.8 s 1Q = 0.9 s n = 1.0. 
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RESULTS AND DISCUSSION 


The results of the sample calculations are presented in figures 5 to 13. In fig- 
ures 5 to 10, the predicted temperature distributions from the present analysis are com- 
pared with exact solutions from the appendix (eqs. (A4) to (A12)). The curves in fig- 
ures 5 to 10 which represent the exact solutions were obtained by use of computer 
programs which evaluate equations (A4) to (A12) at a large number of points. Fig- 
ures 11 to 13 contain only results from the present method since no exact solution was 
available for comparison. 

The agreement between results of the present method and the exact solutions for 
the first four configurations (figs. 5 to 8) is excellent. The numerical values were in 
agreement to eight significant figures; thus, the corresponding curves are coincident. 

For the fifth configuration (fig. 9) using an idealization of 10 equal elements, the 
present results and the exact solution are coincident except in the vicinity of the center 
of the cylinder where the present solution deviates from the exact solution. This slight 
error is attributed to the inability of the present solution to predict the slope discon- 
tinuity exhibited by the exact solution at. the center of the cylinder. It will be recalled 
that one of the assumptions made in the present formulation was that the temperature 
distributions have continuous derivatives. 

In order to show that errors associated with smoothing of slope discontinuities may 
be made negligible, a second calculation was made. In this calculation 20 elements were 
used, smaller elements being placed in the vicinity of the center of the cylinder than at 
the ends. The result of this calculation was a temperature distribution which was coinci- 
dent with the exact solution for plotting purposes. (See fig. 9.) In that figure the 20 ele- 
ment calculation and exact solution are both represented by the solid line. 

The results for configuration 6 are presented in table IV and. figure 10. The results 
were obtained as follows: Using ranges of n of zero to two, zero to three, and zero to 
four, the computer program of block diagram 2 computed the functions u cn (s) and 
u gn (s). All the u gn (s) functions were found to be identically zero and thus agreed with 
the exact solution in equation (A12). The u cn (s) functions were constants and are tab- 
ulated in table IV along with the exact Fourier components from the appendix. 
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Temperature component, 


Angle from strip, radians 


Figure 10.- Temperature distribution in a cylindrical shell with an axial 
strip held at constant temperature. 



0 .2 .4 .6 .8 1.0 

Distance, s 


Figure 11.- Fourier components of temperature distribution in a cylinder 
■with specified temperature at a single point. 




Temperature 




TABLE IV.- FOURIER COSINE HARMONICS OF TEMPERATURE FOR 
CYLINDER WITH GENERATOR HELD AT UNIT TEMPERATURE 


Fourier 

harmonic 

Present 
analysis 
(3 terms) 

Present 
analysis 
(4 terms) 

Present 
analysis 
(5 terms) 

Exact 
solution 
(eq. (A12)) 

n = 0 

0.417 

0.385 

0.368 

0.318 

n = 1 

.417 

.385 

.368 

.318 

n = 2 

.167 

.153 

.147 

.127 

n = 3 

— 

.0769 

.0736 

.0636 

n = 4 

— 


.0433 

.0374 


Inspection of table IV shows that each Fourier harmonic from even the best calcu- 
lation is higher than the corresponding exact harmonic by 16 percent. This error is much 
larger than had been encountered on any of the previous calculations because approxima- 
tions had to be made in this problem which limit the accuracy. These approximations 
are as follows: 

(1) The condition that the entire generator is held at the prescribed temperature is 
approximated by the condition that a finite number of discrete points are held at the pre- 
scribed temperature. 

(2) In contrast to the previous calculations, all the Fourier harmonics must be com- 
puted simultaneously. (See eq. (56).) The storage limitations of the computer require a 
truncation of the number of harmonics which may be computed in a given calculation. 

This truncation necessarily limits the accuracy of the computed harmonics because 
theoretically the exact solution would be achieved only when an infinite number of har- 
monics were computed simultaneously. 

As an additional aid to evaluating the results of this calculation, the Fourier series 
representations of the temperature were computed for two sets of coefficients in table IV, 
the coefficients from the present five-term calculation and the exact coefficients. These 
representations are plotted in figure 10 along with the exact (converged) solution. In 
figure 10 the short-dashed line is the five-term representation based on the present 
analysis. The short-long-dashed line is the five-term representation using exact coeffi- 
cients. The solid line is the converged solution of equation (A10). The hump in the vicin- 
ity of 6 » 7 T in both five-term solutions is characteristic of a truncated Fourier series. 

It can be seen from figure 10 that the result from the present analysis is in reason- 
ably good agreement with the two other results. The present solution is an upper bound 
to both of the other solutions. The maximum disparity between the present five-term 
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solution and the exact five-term representation is 16 percent at 0=0. The maximum 
disparity between the present five-term solution and the converged solution is 29 percent 
at 6 = 7t/4. The maximum difference between the exact five-term representation and 
the converged solution is 16 percent at 0 = 0. 

The final sample calculation is the temperature distribution in a cylinder wherein 
a single point is held at unit temperature. The results are shown in figures 11 to 13. 

This problem was chosen to demonstrate the applicability of the present method to a 
problem for which an exact solution was not available and hence a numerical solution 
was necessary. The point on the cylinder at which temperature is specified is taken to 
be at the midlength of the cylinder at 0=0. The temperature distribution would then 
be expected to be symmetric about the line s = L/2 and the generator 0 = 0. 

The functions u cn (s) and u gn (s) were computed for a range of n of zero to 
four. All the u sn (s) functions were identically zero; thus, the symmetry of the tem- 
perature about 0 = 0 is verified. In addition, as shown in figure 11, the u cn (s) func- 
tions are symmetric about s = L/2. 

In order to illustrate further the nature of the temperature distribution, two pro- 
files through the heated point are calculated. The temperature profile along the generator 
0 = 0 is shown in figure 12. This curve was calculated by adding the curves of figure 11. 
The profile along the circumference s = 0.5 is shown in figure 13. This curve was 
constructed from the equation 

4 

u(O.5,0) = u cn (0.5) cos n0 (57) 

n=0 

The appearance of the curve in the vicinity of 0 = tt/ 2 and 0 = - 17/2 is characteristic 
of a truncated Fourier series. As the number of Fourier terms is increased, the irreg- 
ularities in the curve should disappear. 

The computer programs are very efficient as evidenced by the fact that the total 
running time required on the Control Data 6600 computer system to perform all the cal- 
culations in the paper was less than 2 minutes. 

CONCLUDING REMARKS 

An analytical procedure based on the finite-element method has been developed 
for computing steady-state temperature distributions in thin shells of revolution. The 
shells may have general meridional curvature and nonhomogeneous thermal properties. 
The procedure is applicable to shells subjected to a very general set of thermal loading 
conditions including applied heat flux at the surfaces , convective heat transfer to the 
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surroundings, and equivalent loadings due to specified temperatures on the shell. Gen- 
eral descriptions of two computer programs based on the procedure are given. The first 
program deals with problems wherein temperatures are specified on complete circum- 
ferences only. The second deals with problems wherein temperatures are specified both 
on circumferences and at discrete points. 

The computer programs have been applied to a number of configurations for which 
exact analytical solutions were available: 

(1) An annular plate with specified temperature at the inner radius and insulated 

at the outer radius 

(2) A cylindrical shell with edges at fixed temperatures and with convection to 

surroundings 

(3) A conical frustum subjected to a nonuniform heat flux over the outer surface 

(4) A truncated hemisphere subjected to uniform heat flux at the outer surface 

(5) A cylindrical shell with circumference at midlength held at constant temperature 

(6) A cylindrical shell with a generator held at constant temperature 

In addition to these configurations, a seventh configuration was analyzed for which no 
other solution was available. This configuration was a cylindrical shell in which a single 
point was held at a prescribed temperature. In general, the temperature distributions 
predicted by the present procedure were in excellent agreement with the exact solutions. 
Inaccuracies occurred in the solution for configurations 6 and 7 in which it was theoreti- 
cally necessary to solve for all the Fourier components of temperature simultaneously. 
Thus, a truncation in the size of the set of equations solved was necessary. 

In problems wherein the exact temperature distribution has a slope discontinuity 
such as configuration 5, the present method predicts a continuous slope. In this case 
by using a refined representation of the shell in the neighborhood of the discontinuity, the 
error associated with this phenomenon was made negligible. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., December 28, 1970. 
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APPENDIX 


EXACT SOLUTIONS FOR TEMPERATURE DISTRIBUTIONS 


Taking the first variation of the functional of equation (1) and setting it equal zero 
in the usual manner of the calculus of variations gives the required boundary-value 
problem: 


£ ( K S) + K T § + § ■ 0 “o + H i) u + H o u o + H I U I + <*o + q i = 0 


(Al) 


K|| + h p q F -h F H F (u-u F ) = 0" 
u prescribed 


(s = s 0 ) (A2) 


K S - h L«L + » L H L (» - »l) - 

u prescribed , 


(s = L) (A3) 


Solutions to equations (Al) to (A3) are obtained for the first six shell configurations. 


Annular Disk 

The terms in equation (Al) appropriate to an annular disk with constant thermal 
conductivity subjected to convection to uniform surroundings at zero temperature are 


r = s 


r' = 1 


W 


= 0 


Hq = Hj = H = Constant 


u 0 = uj = 0 
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APPENDIX — Continued 


Qq = q r = 0 


K = Constant 


Equation (Al) becomes 


lii + Idu _ 2 H u = 0 

<j s 2 s ds K 


The boundary conditions are taken to be 


“f = T ° 


du 

ds 


s=l 


= 0 


The solution to the preceding boundary-value problem is given in reference 8 as 

I Q (Xs) Kj(X) + K 0 (XS) I 1 (X) 


where 


*1 


x 2 = 2 h 

K 


u(s) = T 0 


[ o(|) K iW + K o(|)iiW 


modified Bessel function of the first kind of order zero 


modified Bessel function of the first kind of order one 


modified Bessel function of the second kind of order zero 


modified Bessel function of the second kind of order one 


(A4) 


Numerical results were obtained for the following values of the parameters: 


K = T 0 = 1 
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APPENDIX - Continued 

Cylindrical Shell With Convection to Constant Ambient Temperature 

For the cylindrical shell with convection to constant ambient temperature, when 
constant properties of the shell and surroundings are assumed, the terms in equa- 
tion (Al) are 

K = Constant 

r = Constant 

% = Constant 

Hj = 0 

QO = qj = 0 

A =0 

de 

Uq = Constant 


The governing differential equation is reduced to 


H 2 

K — - H n u + H n u n = 0 
ds 2 ° ° ° 


The boundary conditions are taken to be 
u(0) = u(L) = 0 
The solution is 


u(s) = u^fl - cosh Xs) + u J cosh ^ ~ A sinh Xs 
° v J °\ sinh XL J 

where 

x 2 = ! 1 ° 

K 


Equation (A5) was evaluated for the following values of the parameters: 


u o = H o = K = L = 1 


(A5) 


Conical Frustum Subjected to Nonaxisymmetric Heat Flux at Outer Surface 

For the conical frustum subjected to nonaxisymmetric heat flux at the outer surface, 
the terms in equations (Al) to (A3) are 
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APPENDIX - Continued 


K = Constant 
r(s) = -j^s 

2 

Ho = Hj = 0 

qo = Q(0) 
q x = o 

u(s 0 ,e) = T 0 (e) 

u(L,e) = T L (fi) 

Equations (Al) to (A3) are 

9 2 u , l.ga , _4__ 9 2 u , Q(8) _ Q 
8s 2 s 8s 3s 2 de 2 K 

u(s 0 ,e) = T o (0) 
u(L,e) = t l (s) 

Consider the case: 

Q(0) = Qn cos n0 
T o (0) = T Q n cos nfl 

T L (0) = T L?n cos ne 

It follows that 

u(s,0) = u n (s) cos n 6 

thus, 

2 

d u n I du n 4n 2 Q n _ n 
ds 2 s ds ' 3s 2 n K 

u n( s o) = T 0?n 
u n( L ) = T L,n 
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APPENDIX - Continued 


The solution for n =£ 0 is 


where 


u n( s ) = A n s Mn + B n s ^ 


Qn s 

k(m„ 2 - 4) 


(A6) 


f x n = — n 

V 3 


An ~ 


L" Mn 

T - Q ” S ° 2 " 

- s " Mn 

T V* 2 


_ ,IX K( Mn 2 -4) 

s o 

_ 1 ’ n K(/x n 2 - 4)_ 


Mn T ~^n a “Mn T ^n 
s o L " s o L 


B n = 


-L _Mn 

T _ Q n s o 




Q n L 2 " 


L u "‘ k (»*„ 2 - *i 

i, n K) 

k 2 - 4 


s /n L -^n _ ^n L ^n 


For n = 0, equations (Al) to (A3) become 

d 2u O 1 du Q Q q 
ds 2 s ds K 

u 0 ( s o) = T o 
u(L) = T L 

The solution is 


u G (s) = A 0 + B 0 log e s - 


4K 


where 


Q 0 s o 


A 0 = 


ioge L\T 0 + - log e s 0 ^T l + 


Qq^ 

4K 


l°gp (L/s 0 


(A7) 
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APPENDIX - Continued 


Q 0 L <W 

T t H — — “ Tq °° 


B 0 = 


4K 


4K 


log e( L / s o) 

Numerical results were obtained for the following values of the parameters: 


s o = 1 


L = 10 

Qn=l 


T 0,n = 1 


Ti n = 6 

l,n 


K = 1 


(n = 0, 1, . . 5) 


Truncated Sphere Subjected to Axisymmetric Heat Flux at Outer Surface 

For the truncated sphere subjected to axisymmetric heat flux at the outer surface, 
the appropriate terms in equations (Al) to (A3) are 

K = Constant 

r = R cos(s/R) 

r' = -sin(s/R) 

qj = 0 

q^ = Q = Constant 



u(0) = T 0 
u(L) = T l 

where R is the radius of the sphere. Equation (Al) takes the form 
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APPENDIX — Continued 


drj 2 K 


where 


77 = sin — 
R 


The general solution with 77 as the independent variable is 
u(t 7 ) = £ log + B + log(l - ? 7 2 ) 

Z 1 - 7} 

The solution with s as the independent variable is 

u(s) = A l0 Jl ±jjn(gM| + B ♦ 95! logfcosfljl 
2 [l - sin(s/R)J K L VR7j 


where 


2<T L -T 0 -£3-log 


A = 


log 


1 + sin(L/R) 
1 - sin(L/R) 


[ cos (t) 

i 

)J 


B = T f 


(A8) 


Numerical results were obtained by evaluating equation (A8) for the following values 
of the parameters: 

R = 1 

L = 1 



Q = 1 


Cylinder With Circumference at Midlength Held at Fixed Temperature 

The solution for the cylinder with circumference at midlength held at fixed tempera- 
ture is given in reference 9 as 
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APPENDIX — Continued 


V 0" V 1 

Vi + 77 rr\ sinh ix s 


1 ^W) 


u(s) = ( 


V 0 - Vj 

Vi + — u a.tv Sinh ju(L - s) 


1 ^ 




(A9) 


where 



Vj temperature of edges of cylinder 

Vq temperature of circumference at midspan of cylinder 

L length of cylinder 

Numerical results were obtained by evaluating equation (A9) for the following values of 
the parameters: 


Ho = K=l 
Vi = i 

v 0 = io 

L = 1 


Cylinder With a Generator Held at Fixed Temperature 

The exact solution for the cylinder with a generator held at fixed temperature is 
given in reference 9. The solution is: 


where 


u(e) = t cosh M<g - Hi 

COSh JU.7T 


(A10) 




Ho 

K 


T temperature of generator 
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APPENDIX — Concluded 


For comparing this solution with results of the present analysis, the Fourier series 
components of the expression in equation (A10) are computed. 

Thus, 

N 

u(0) = ^ cos nSH-^bn sin n6 (All) 

n=0 

where 


T 

a A = — tanh ju.7r 

0 l±TI 


a n 


2T 


77 (l + n 2 ) 


tanh ju.7r 




b n = 0 




(A12) 


For the numerical calculations, 


T = Hq = K = 1 

thus, 

M = 1 


L-7324 
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